/*<gaussnode.cu
 * GPU-accelerated test program to determine nodes and weights.
 * In future will be included as main method for nodes and weight creation.
 * Copyright (C) 2011 Ivan Medvedev
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
#include <stdio.h>
#include <cuda.h>
#include <math.h>
#include <time.h>
#define ACC 16
__global__ void
gaussnodes(float* nodes,float* net, float* weights,int N, float start, float end) {
float step;
float wg[16] = {0.027152459411754094852, 0.062253523938647892863, 0.095158511682492784810, 0.124628971255533872052, 0.149595988816576732081, 0.169156519395002538189, 0.182603415044923588867, 0.189450610455068496285, 0.189450610455068496285, 0.182603415044923588867, 0.169156519395002538189, 0.149595988816576732081, 0.124628971255533872052, 0.095158511682492784810, 0.062253523938647892863, 0.027152459411754094852  };
float xg[16] = { 0.989400934991649923596, 0.944575023073232576078, 0.865631202387831743880, 0.755404408355003033895, 0.617876244402643748447, 0.458016777657227386342, 0.281603550779258913230, 0.095012509837637440185, 0.095012509837637440185, 0.281603550779258913230, 0.458016777657227386342, 0.617876244402643748447, 0.755404408355003033895, 0.865631202387831743880, 0.944575023073232576078, 0.98940093499164992359  };
//__shared__ float net[4];
step = (end-start)/(float)N;
/*int shift;
for (int vec=blockIdx.x; vec< N/16+1; vec+=gridDim.x){
	for(int ng=threadIdx.x; ng<ACC; ng+=blockDim.x){
	weights[ng+vec*16]=ng;
	}
}
*/
//run on every pair of input vectors; each pair is in a different block in the grid
for (int vec_n=blockIdx.x; vec_n < N; vec_n += gridDim.x)
{
    //run on every slice of input vectors
    for (int slice =threadIdx.x; slice < ACC; slice += blockDim.x)
    {
        //each thread calculates one product
        weights[slice+vec_n*gridDim.x] = step*wg[slice]*0.5;
        __syncthreads();
}
}
int i = blockDim.x * blockIdx.x + threadIdx.x;
	if(i<N+1){
	net[i]=start+i*step;
	}
for (int vec_n=blockIdx.x; vec_n < N; vec_n += gridDim.x)
{
//nodes[vec_n] = 778;
    //run on every slice of input vectors
    for (int slice =threadIdx.x; slice < ACC; slice += blockDim.x)
    {
	float mgn=net[vec_n]+net[vec_n+1];
        //each thread calculates one product
	if(slice/2==slice/2.){
        nodes[vec_n*gridDim.x+slice] = 0.5*(mgn+step*xg[slice]);
	}else{
	nodes[vec_n*gridDim.x+slice] =0.5*(mgn-step*xg[slice-1]);
	}

        __syncthreads();
}
}

}

int main(int argc, char** argv)
{

float* h_X;
float* d_X;
float* h_Y;
float* d_Y;
float* h_Z;
float* d_Z;
float start,end;
start=0;
end=10;
h_X=NULL;
d_X=NULL;
h_Y=NULL;
d_Y=NULL;
int N=10; //size of vector
    // Allocate vectors in device memory
size_t size = N*16*sizeof(float);
h_X = (float*)malloc(size);
h_Y = (float*)malloc(size);
cudaMallocHost( (void**)&h_X, size);
cudaMalloc((void**)&d_X, size);
cudaMallocHost( (void**)&h_Y, size);
cudaMalloc((void**)&d_Y, size);
h_Z = (float*)malloc(size/16+1);
cudaMallocHost( (void**)&h_Z, size);
cudaMalloc((void**)&d_Z, size);
    // Copy vectors from host memory to device memory
cudaMemcpy(d_X, h_X, size, cudaMemcpyHostToDevice) ;
    // Invoke kernel
//cudaMemcpyAsync(d_X, h_X, size,cudaMemcpyDeviceToHost, 0) ;
    int threadsPerBlock =16;
    int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;
cudaThreadSynchronize();
    gaussnodes<<<threadsPerBlock,blocksPerGrid>>>(d_Y,d_Z,d_X,N,start,end);
//    cutilCheckMsg("kernel launch failure");
cudaMemcpy(h_X, d_X, size, cudaMemcpyDeviceToHost) ;
cudaMemcpy(h_Y, d_Y, size, cudaMemcpyDeviceToHost) ;
for(int i=0; i!=20; i++){
	printf("%f\n",h_Y[i]);
}
cudaFree(h_X);
cudaFree(d_X);
cudaFree(h_Y);
cudaFree(d_Y);
h_X=NULL;
d_X=NULL;
return 0;

}	
